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For strongly repulsive bosons in one dimension, we provide detailed modeling of the Bragg pulse 
used in quantum Newton’s cradle-like settings or in Bragg spectroscopy experiments. By employing 
the Fermi-Bose mapping for a finite harmonically trapped gas and the Quench Action approach for 
a thermodynamic system on a ring, we reconstruct the exact post-pulse many-body time evolution 
of Lieb-Liniger gases in the Tonks-Girardeau limit, together with their changing local density profile 
and momentum distribution. Our results display a clear separation of timescales between rapid and 
trap-insensitive relaxation immediately after the pulse, followed by slow in-trap periodic behaviour. 


The study of many-body quantum physics has in re¬ 
cent years been transformed by the progress achieved in 
experiments on ultracold atoms [T]. The context of one¬ 
dimensional (ID) bosonic gases in particular provides a 
fertile ground for investigating physics beyond traditional 
paradigms [2] , with concepts such as the Luttinger liquid 
and exact solvability [3] playing a primary role. 

One of the main experimental probes of cold gases is 
Bragg spectroscopy HE], which consists in applying a 
pulsed monochromatic laser grating onto the gas, thereby 
creating excitations at (multiples of) the recoil momen¬ 
tum q, where q/2 is the wavevector of the laser. The 
precise time-dependent form of the pulse can be chosen 
to optimize (de)population of specific quantum states. 
In [II [8] , a two-pulse sequence was used to optimize the 
population of the first Xq satellites as compared to the 
zero-momentum ground state of a Bose-Einstein conden¬ 
sate. The theoretical description of such a sequence relied 
on a two-state model in which many-body dynamics are 
not included. In one dimension however, many-body ef¬ 
fects are inescapable. One of the fundamental models in 
this context is the Lieb-Liniger gas [9] of bosons in a ID 
continuum interacting with a ^-function potential, giv¬ 
ing a proper description of atoms in tight transverse con¬ 
finement m- This model is relevant to the description 
of experiments, most prominently the famous quantum 
Newton’s cradle experiment m, in which a Bragg pulse 
is used to initiate the oscillations. Bragg spectroscopy 
has also recently been used to investigate correlated ID 
Bose gases of rubidium [12] and cesium m- In these, 
the heating of the gas resulting from the Bragg pulse was 
measured and matched using linear response to theoret¬ 
ical calculations of the dynamical structure factor of the 
Lieb-Liniger gas at finite temperature El- 

Our main objective here is to model the effects of Bragg 
pulses theoretically for correlated ID Bose gases, from 
first principles and without approximation (and thus be¬ 
yond linear response), for experimentally relevant se¬ 
tups. We focus on the Tonks-Girardeau limit UMlj of 
strongly repulsive Lieb-Liniger bosons both on a periodic 
interval and in confining traps [T8ll26] . Instantaneous 


Bragg pulses of varying amplitude A and wavevector q 
are studied via their effect on physical observables: the 
time-dependent local density of the gas, and the exper¬ 
imentally more easily accessible momentum distribution 
function (MDE). 

We start by modeling the Bragg pulse as a standing 
wave forming a one-body potential V{x) = Vb cos(gx) 
that couples to the density p{x) = T’f'(x)T(x), where the 
Bose fields obey the canonical equal-time commutation 
relations, [T(x), T’f'(^)] = S{x — y). Eor a general Bragg 
pulse the gas is perturbed for a finite duration Tq. We 
will however consider the regime where the motion of the 
particles during the pulse can be neglected (the Raman- 
Nath limit), in which case the Bragg pulse is also referred 
to as a Kapitza-Dirac pulse EZHUj. Taking the limit 
To ^ 0 such that A = VqTq is kept finite, the Bragg pulse 
operator Ub given by 

= , 1 ) 

where we have used the convention h = 1. Applying 
such an instantaneous Bragg pulse to a ground state 
IV’Gs) yields a post-pulse state I'tpq^A) = 17B(g, ^)|V’Gs), 
which can be interpreted as the initial state of a quantum 
quench [29H3T] . Typical experimental pulses HMUsg 
correspond to Bragg momentum q ~ 27rn and A ~ 1, 
where n is the mean density. 

The post-pulse time evolution is driven by the Lieb- 
Liniger (LL) model of interacting bosons 

AT ^2 

= T. ( 2 ) 

i=l * l<i<j<N 

In what follows we will focus on the hard-core Tonks- 
Girardeau (TG) limit c ^ oc [T5j|T6] and we will consider 
the model both on a ring (periodic boundary conditions) 
and in a parabolic trapping potential. 

In the hard-core limit, the bosonic many-body 
wavefunction can be related through the Eermi-Bose 
(EB) mapping [16] to the many-body wavefunction 


2 


of free fermions = ni<i<j<Ar sgn(xi - 

Xj) 2 pF{x;t), where x and the fermionic 

wavefunction is the usual Slater determinant of the 
free single-particle (SP) wavefunctions, = 

detw [tl^j{xi]t)] /y/N~\. Following [T9l|24l[2^, the bosonic 
one-body density matrix, defined as t)) = 

A^/dx 2 ...dxAr is given in 

terms of a single determinant involving the time- 
dependent fermionic SP states. This allows for 
an efficient computation of the MDF (h(/c,t)) = 

^ Jdxdy 

Starting with the ring geometry, our ground state con¬ 
sists of SP plane waves, on which the Bragg pulse im¬ 
prints a cosine phase due to the one-body potential, 

0) = (3) 

vL 

with ground-state rapidities ^ j)}f=i 

forming a Fermi sea with Fermi momentum Xp = 

Note that the Bragg momentum is quantized due to the 
periodic boundary conditions of the ring, q = ^riq with 
Uq e N. Expanding Eq. in plane waves, the time- 
dependent SP wavefunctions after the Bragg pulse yield 


for all Of G N. This leads to the following stationary-state 
distribution [31 EH, 

p 7A^) = [0{X-l3q+XF)-0{X-/3q-XF)]MiA)f 

( 7 ) 

where 0 is the Heaviside step function. The saddle point 
distribution is a sum of copies of the ground-state density 
of rapidities, pgs(A) = ^ [^(A +Ai?) — ^(A —Ai?)], shifted 
by multiples of q and weighted by the modified Bessel 
functions. 

This form of the stationary state is consistent with the 
QA approach [36] , which furthermore states that the time 
evolution of local observables after a quantum quench 
is given by a sum over particle-hole excitations in the 
vicinity of Pg^^(A) |34l |38]. One easily obtains the 

time evolution of the density of the gas, 

T77T7 

limth {i>q,A (t) I \'lpg,A (t)) = X (8) 

^ Jp{-2Asm{q'^/3t/2m)) , 

/3 = —CX3 


'tpj{x-,t) 


oo ^ 

Ip(^—iA) _ ^-i{Xj+/3q)x^-i{Xj+/3q)^t/2m ^ 

( 4 ) 


where Ip{z) is the modified Bessel function of the first 
kind. 

Contrary to the finite-size Eermi-Bose mapping, the 
Generalized Gibbs Ensemble (GGE) |30l |33] and the 
Quench Action (QA) approach [34l [35] enable the study 
of the Bragg pulsed system (on a ring) in the thermody¬ 
namic limit {N oc with N/L fixed). The GGE can 
be constructed using the infinite number of conserved 
charges {Qa}^=i provided by the integrability of the 
LL model, with Q 2 = 2mH^ and eigenvalues Qa(A) = 
associated to a Bethe state |A) = |Ai,..., Aat). 
The expectation values of the charges on the initial post¬ 
pulse state can be computed using the overlaps {Xl^pq^A)^ 
which can be derived from the matrix elements for the 
Bragg pulse between two Bethe states |A) and \/x) [36] 


{fx\UB{q^A)\X) 


= det^v 




Xj 


( 5 ) 


where we defined 5^7 = 5(x-^l,)mod qfi- Taking the 
thermodynamic limit, the energy density of the system 
after the Bragg pulse yields limth ^ {^qA \ Q 2 = 

|A|. + \ {qAY. The GGE logic |30l|33] then requires the 
expectation values of all charges to be reproduced by the 
equilibrated post-pulse system, described by a density of 
rapidites Pg^y[(A), i.e. 

limth J (V’q.Al Qa Hq,A) = / <^X A“ , (6) 

J —OO 


with Jj 3 {z) the Bessel function of the first kind. The time 
evolution of the density is compared to the EB result 
for A" = 50 in Eig. and shows excellent agreement, 
with relative differences of order 0.4% due to finite-size 
effects in the EB calculations. Throughout the paper, all 
data is produced setting m = 1. As a consequence of 
the Raman-Nath limit, the initial post-pulse density (at 
t = 0) is unaltered from the fiat ground state profile. A 
sharp density profile then develops, mimicking the one- 
body cosine potential that was instantaneously turned 
on and off, followed by relaxation back to a flat profile at 
time scales t ~ m/qXp = {qvs)~^, with Vg = Trnjm the 
sound velocity. 



FIG. 1. Time evolution of the density after a Bragg pulse 
with q — and A = 1.5, computed by (a) the FB mapping 
and (b) the QA approach. The relative differences between 
the two results due to finite-size effects are less than 0.4%. 
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The QA approach also provides access to the time evo¬ 
lution of the MDF [36l|39]. The result is plotted in Fig.|^ 
along with the FB result for N = 50. Except for mi¬ 
nor disagreements in the sharp peaks due to finite-size 
effects, the large-system-size dynamics after the Bragg 
pulse is again well captured by a Fermi-Bose mapping 
for = 50 particles. At t = 0, using the commutation 
relations between the Bragg operator and the Bose fields, 
one can show that the momentum distribution is simply 
a sum of copies of the ground-state MDF [36], with a 
small-/c divergence {n{k))GS ^ centered around 

multiples of q. The sharply peaked MDF then relaxes to 
a characteristic ghost-like shape m, with the mixing of 
particles with different momenta causing a substantially 
broadened stationary MDF. 

In Fig. the equilibrated MDF is shown for differ¬ 
ent values of q and A. Similar to the initial MDF, the 
late-time distribution behaves like a superposition of in¬ 
dependent peaks shifted to multiples of q. The width of 
each satellite shows no dependence on the value of q, and 
is only influenced by the choice of A. Since in the limit of 
A ^ 0 the resulting MDF is just that of the ground state 
and would stay constant as time progresses, the broad¬ 
ening can be ascribed to interactions between particles 
belonging to different copies of the ground state density 
of rapidities in 



FIG. 2. Time evolution of the MDF after a Bragg pulse with 
q — and A = 1.4, computed with the QA approach (left 
half) and the FB mapping (right half). Because the FB map¬ 
ping treats a finite system (N = 50) the momenta are quan¬ 
tized, causing less pronounced peaks for short times. All other 
results are in excellent agreement with the QA computations. 

Next, we will use the FB mapping to investigate how 
these observations translate to the more experimentally 
relevant geometry of a harmonic trapping potential, with 
the Hamiltonian i^trap = ^ll + and cj 

the trapping frequency. The ground state SP harmonic 
oscillator wavefunctions are given by 

'4’jix) = -j= (^ —j e 2 Hj {^/mxx) , (9) 



FIG. 3. The relaxed MDF function for N = 50 particles as 
a function of the Bragg momentum q. In the upper panel 
A = 1, and in the lower panel A is fixed to 1.5. 

for j = 1,..., A", with Hj{x) denoting the Hermite poly¬ 
nomials. Similar to Eq. ([^, acting on these states with 
the Bragg operator Ub leads to an additional phase of the 
one-body cosine potential in the SP wavefunctions. Us¬ 
ing the propagator for the quantum harmonic oscillator 
(the Mehler kernel) |40|, we compute the time evolution 
of the SP wavefunctions [36] : 

oo 

/3= —OO 

( 10 ) 

As a consequence of the Tonks-Girardeau limit, the SP 
wavefunctions are periodic in time with period 27r/ci;, 
which is reflected in observables such as the density and 
the MDE. This periodicity is expected to be broken by 
finite-c interactions and anharmonicities in the trapping 
potential. The time evolution of the density and the 
MDE during one period is shown in Eig. [^ where the con¬ 
tributions from particles belonging to different momen¬ 
tum satellites are clearly distinguishable. During the ini¬ 
tial stages of relaxation (and around multiples of t = ^) 
the density shows strong oscillations and the initially 
sharply peaked MDE relaxes rapidly to a more broad¬ 
ened shape. This prerelaxation is well separated from 
the collective periodic motion due to the trap, suggest¬ 
ing that it is governed by the same physics as relaxation 
on a ring. 

In Eig. [^the density at early stages in the oscillation 
cycle is compared with the density dynamics on a ring, 
where the latter was supplemented by a local density ap¬ 
proximation (EDA) to account for the classical expansion 
of the gas in the trap [ 31111 ]. The initial density profile 
is accurately reproduced by the LDA, except for small 
differences near the edges originating from the gradient 
in the local density not accounted for by the LDA [42ll45] . 
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FIG. 4. The time evolution of the density (top) and MDF 
(bottom) in the trap, computed with the FB mapping for 
iV = 50, a; = 10/iV, A = 1.5 and q = Stt. 


Note however that these differences do not stay confined 
to the edges and propagate towards the center as time 
progresses. 

The short-time MDF in the trap and ring geometry 
is shown in Fig. up to t = 0.02257r/cc; . The initial 
distributions are nearly identical, after which the MDFs 
dephase in a similar fashion to a (pre)relaxed ghost-like 
shape. The strong similarities can be attributed to the 
short-range correlations characterizing the post-quench 
steady state, for which the one-body density matrix de¬ 
cays exponentially as a function of the distance between 
the particles. Large-distance effects due to the trap ge¬ 
ometry lead to discrepancies only at low momenta k. We 
conclude that the short-time dynamics in a trap closely 
resembles the dephasing on a ring and is thus governed by 
the physics of hard-core interactions. The time scale as¬ 
sociated to this (pre)relaxation is much shorter than the 
collective oscillations in the trap. Considering conditions 
similar to the Newton’s cradle experiment, we estimate 
the short time scale to be of the order of 10 /rs. Since this 
estimate is of the same order of magnitude as the pulse 
durations used in m, an interesting next step would be 
to extend the study of the effects of interactions to longer 
pulses in the Bragg regime. 



FIG. 5. Time evolution of the density in a trap, computed 
with: the FB mapping for N = 50 particles (left), the QA ap¬ 
proach on a ring with an LDA accounting for the trap (mid¬ 
dle). The difference between the two results is shown on the 
right. The Bragg pulse parameters are set to A = 1.5 and 
g = TT with a trapping frequency u = 10/N. 
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FIG. 6. Time evolution of the MDF for the trap geometry 
(solid lines) and the ring geometry (dashed lines), obtained 
with the FB mapping for A" = 50 particles. The trapping 
frequency is set to uj = 10/A, and the Bragg parameters are 
given by A = 1.5 and q = Stt. The time step At is set to 


An interesting open question is how these observations 
extend to finite-c interactions. Away from the TG limit, 
the initial post-pulse MDF remains a weighted sum of 
copies of the ground-state MDF [36] , leading to a decreas¬ 
ing spread of momenta in the satellites centered around 
multiples of q, as one goes from the hardcore limit to the 
BEC limit (c ^ 0). Although finite-c dynamics is cur¬ 
rently unattainable, we expect the time scale of the rapid 
relaxation to grow for smaller c, proportional to the in¬ 
verse of the sound velocity, Vs = ^ A- + ...). 

Conclusion. In summary, we have developed a theo¬ 
retical description of the Bragg pulse for one-dimensional 
Bose gases and shown that the time evolution of physi¬ 
cal observables for a Bragg pulsed Tonks-Girardeau gas 
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in a trap is characterized by two well-separated time 
scales. The shortest time scale is dominated by the trap- 
insensitive hard-core interactions and causes a substan¬ 
tial broadening of the momentum distribution function 
well before the collective motion due to the presence of 
the trap sets in. Our work opens up the possibility to 
study the influence of interactions on more general pulse 
protocols and their detailed effects on experimentally rel¬ 
evant observables. 

We thank M. Brockmann, N.J. van Druten, V. Grit- 
sev, F. Meinert, H.-C. Nagerl, J. Schmiedmayer, F.E. 
Schreck, D. Weiss and J. van Wezel for useful discussions. 
This work was supported by the Netherlands Organisa¬ 
tion for Scientific Research (NWO) and the Foundation 
for Fundamental Research on Matter (FOM), and forms 
part of the activities of the Delta-Institute for Theoreti¬ 
cal Physics (D-ITP). This research was done in part un¬ 
der the auspices of the CMPMS Dept, at Brookhaven 
National Laboratory, which in turn is supported by the 
U.S. Department of Energy, Office of Basic Energy Sci¬ 
ences, under Contract No. DE-AC02-98CH10886. We 
are grateful for support from the Centre de Recherches 
Mathematiques of the U. de Montreal, where this work 
was completed. 


* R.vandenBerg2@uva.nl 



[3] T. Giamarchi, Quantum Physics in One Dimension (Ox¬ 
ford University Press, 2004). 


[4] P. J. Martin, B. G. Oldaker, A. H. Miklich, and D. E. 
Pritchard, Phys. Rev. Lett. 60, 515 (1988) 

[5] J. Stenger, S. Inouye, A. P. Chikkatur, D. M. Stamper- 
Kurn, D. E. Pritchard, and W. Ketterle, Phys. Rev. 
Lett. 82, 4569 (1999) 

[6] R. Ozeri, N. Katz, J. Steinhauer, and N. Davidson, Rev. 
Mod. Phys. 77, 187 (2005) 

[7] Y.-J. Wang, D. Z. Anderson, V. M. Bright, E. A. Cornell, 
Q. Diot, T. Kishimoto, M. Prentiss, R. A. Saravanan, 
S. R. Segal, and S. Wu, Phys. Rev. Lett. 94, 090405 
(2005) 

[8] S. Wu, Y.-J. Wang, Q. Diot, and M. Prentiss, Phys. Rev. 
A 71, 043602 (2005) 

[9] E. H. Lieb and W. Liniger, Phys. Rev. 130, 1605 (1963). 

[10] M. Olshanii, Phys. Rev. Lett. 81, 938 (1998) 

[11] T. Kinoshita, T. Wenger, and D. S. Weiss, Nature 440, 
900 (2006) 

[12] N. Fabbri, M. Panfil, D. Clement, L. Fallani, M. Inguscio, 
C. Fort, and J.-S. Caux, Phys. Rev. A 91, 043617 (2015), 

[13] F. Meinert, M. Panfil, M. J. Mark, K. Lauber, J.-S. Caux, 
and H.-C. Nagerl, Phys. Rev. Lett. 115, 085301 (2015) 

[14] M. Panfil and J.-S. Caux, Phys. Rev. A 89, 033605 
(2014) 

[15] L. Tonks, Phys. Rev. 50, 955 (1936) 

[16] M. Girardeau, J. Math. Phys. 1, 516 (1960) 


[17] E. Haller, M. Gustavsson, M. J. Mark, J. G. Danzl, 
R. Hart, G. Pupillo, and H.-C. Nagerl, Science 325, 
1224 (2009) 

[18] A. Minguzzi and D. M. Gangardt, Phys. Rev. Lett. 94, 
240404 (2005) 

[19] R. Pezer and H. Buljan, Phys. Rev. Lett. 98, 240403 
(2007) 

[20] D. M. Gangardt and M. Pustilnik, Phys. Rev. A 77, 
041604 (2008) 

[21] D. Muth, B. Schmidt, and M. Fleischhauer, New Journal 
of Physics 12 , 083065 (2010) 

[22] C. Schenke, A. Minguzzi, and F. W. J. Hekking, Phys. 
Rev. A 84, 053636 (2011) 

[23] G. E. Astrakharchik and L. P. Pitaevskii, Europhys. Lett. 
102, 30004 (2013) 

[24] M. Collura, S. Sotiriadis, and P. Calabrese, Phys. Rev. 
Lett. 110, 245301 (2013) 

[25] E. Quinn and M. Haque, Phys. Rev. A 90, 053609 (2014) 

[26] F. Cartarius, E. Kawasaki, and A. Minguzzi, 
arXiv:1505.01009 

[27] P. L. Kapitza and P. A. M. Dirac, Proceedings of the 
Cambridge Philosophical Society 29, 297 (1933). 

[28] D. L. Freimund, K. Aflatooni, and H. Batelaan, Nature 
413, 142 (2001) 

[29] P. Calabrese and J. Cardy, Phys. Rev. Lett. 96, 136801 
(2006). 

[30] M. Rigol, V. Dunjko, V. Yurovsky, and M. Olshanii, 
Phys. Rev. Lett. 98, 050405 (2007). 

[31] A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalat- 
tore. Rev. Mod. Phys. 83, 863 (2011). 

[32] R. E. Sapiro, R. Zhang, and G. Raithel, Phys. Rev. A 
79, 043630 (2009) 

[33] M. Rigol, V. Dunjko, and M. Olshanii, Nature 452, 854 
(2008). 

[34] J.-S. Caux and F. H. L. Essler, Phys. Rev. Lett. 110, 
257203 (2013). 

[35] J. De Nardis, B. Wouters, M. Brockmann, and J.-S. 
Caux, Phys. Rev. A 89, 033601 (2014). 

[36] See Supplemental Material. 

[37] G. Bevilacqua, V. Biancalana, Y. Dancheva, T. Mansour, 
and L. Moi, J. Math. Phys. 52, 033508 (2011) 

[38] J. De Nardis, L. Piroli, and J.-S. Caux, 
arXiv:1505.03080. 

[39] J. De Nardis and J.-S. Caux, J. Stat. Mech.: Th. Exp. 
2014, P12012 (2014) 

[40] F. G. Mehler, J. Reine Angew. Math. 66, 161 (1866) 

[41] A. S. Campbell, D. M. Gangardt, and K. V. 
Kheruntsyan, Phys. Rev. Lett. 114, 125302 (2015) 

[42] K. V. Kheruntsyan, D. M. Gangardt, P. D. Drummond, 
and G. V. Shlyapnikov, Phys. Rev. A 71, 053615 (2005) 

[43] A. H. van Amerongen, J. J. P. van Es, P. Wicke, K. V. 
Kheruntsyan, and N. J. van Druten, Phys. Rev. Lett. 
100, 090402 (2008) 

[44] J. Armijo, T. Jacqmin, K. Kheruntsyan, and 1. Bou- 
choule, Phys. Rev. A 83, 021605 (2011) 

[45] T. Jacqmin, J. Armijo, T. Berrada, K. V. Kheruntsyan, 
and 1. Bouchoule, Phys. Rev. Lett. 106, 230405 (2011) 







6 


SUPPLEMENTAL MATERIAL 


Matrix elements and initial state on the ring 

The initial state after the Bragg pulse is easily obtained 
from the matrix elements of the operator 

UB{q,A) = I dx cosiqx)if^ {x)if{x) ^ 

In the Tonks-Girardeau (TG) limit on a ring geometry 
we use the eigenstates 

|A) = ^ d^x^pNix\X) ... ^^(xn) |0) , 

( 12 ) 

with wavefunctions given by 
xPn{x\X) = —h det j JJ sgn{xj - xi). (13) 

l<i<j<iV 

By commuting Usiq^A) through the creation operators, 
the matrix elements can be expressed as 

{X\UB{q,A)\^l) = ^ J d^xd^x'i;N{x\X)*'tpN{x'\^l) 

X (0| II 4'(x;) 4-^(xj)|0). 

n 7 

(14) 

The expectation value of the bosonic operators conspires 
with the signs in the Tonks-Girardeau wavefunctions 
leading to a determinant of (5-functions. Treating the 
coordinates as dummy variables under the integral sign, 
it is easy to rewrite the integral in factorized form as a 
determinant of integrals of the form 

(15) 

where we define = (5(a-^) modg,o where we used 

that Xj — lik and q he on the momentum lattice (^'KajV) 
with (a G Z. This results in the matrix elements 

\ « / j,k 

(16) 

The initial state 

\i^g,A)=UB{q,A)\ijGs) (17) 

is easily expressed in the Tonks-Girardeau eigenbasis us¬ 
ing the matrix elements {/j.\UB{q^ A)\2pGs)• 


{fi\UB{A)\X) 

- 7m -= detjv 


The stationary state on a ring from a GGE and the 
Quench Action approach 

In order to implement the GGE logic [30j [33], one 
starts with computing the conserved charges on the ini¬ 
tial state. Let us focus on the case q > 2 Xf, for which 


the overlaps {X\^pq^A) coming from Eq. (16) reduce to a 
simple product of N modified Bessel functions. While 
odd charges are trivially zero, for the even charges we 
find at finite system size 

102a 


N 


V 2a 


= J2J2dpdA)f{xM) 

j=i /3ez 

N a 

j=l 

N a 


1=0 


21 


\qPY 


j=l1=0 


= ( 18 ) 


where we defined Xj{/3) = + qp and where the co¬ 

efficients B2i^o come from the sum over the order of the 
Bessel functions and are known recursively [37] . The sum 
over particles j can be performed, after which the ther¬ 
modynamic limit can be taken. 


li^th lO2a/-A^|'0q,A) 

^ V- Ga\ (A) 

E G?/ 2(a-0 + 1 ’ 


(19) 


where n is the average particle density. Eor example, the 
energy density pumped into the system by an instanta¬ 
neous Bragg pulse is given by 


Ya‘^ 


( 20 ) 


limth {{i’q,A\Q2/N\i;g^A) - (V’gs|Q2/1V|V’gs)) = 

One can show that the saddle-point density 

p 7A^) = [O{X-/3q+XF)-O{X-/3q-XF)]\l0{iA)f 

^ /3ez 

( 21 ) 

reproduces these values of the conserved charges, i.e. 
L I-^^^p7,aW = limth(V’q,A|Q 2 a|'*/’ 9 ,yi) for all a € 

N, by performing the integral and recasting the infinite 
sum into the coefficients 52a,o- Since the local conserved 
charges (if well defined) uniquely determine the saddle 
point, we have thus found the saddle-point density after 
a Bragg pulse for q > 2 Xf. Eor smaller Bragg momenta 
q <2 Xf the computation becomes considerably more dif¬ 
ficult due to the determinant structure of the overlaps, 
but one can show that the saddle-point density given in 
Eq. ([^ is still correct. 

reproduces 


The Quench Action (QA) approach 
this saddle-point density for g' > 2Ai?. As a conse¬ 
quence of working in the Tonks-Girardeau regime, there 
are many microstates with exactly the same overlap. We 
can rephrase the overlaps as 


({A,(/?i)}f=i|^,,A) = L^ n [Ia{-iA)Y 


( 22 ) 
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where is the number of rapidities j with Pj = a and 
a G Z. In the thermodynamic limit these numbers are 
given by 


Lagrange multiplier h is added to fix the particle density 
to n = N/L. The resulting GTBA equation is not an 
integral equation because of the Tonks-Girardeau limit, 


paq-\-XF 

na = L dX p{X) . 

Jaq—Xp 


(23) 


The normalized overlap coefficients = 

— \n(^{{Xj{Pj)}^^i\'ipq^A)/L^) have a well-defined 

thermodynamic limit, 

S[p] = limth Re (24) 

^ POiq-\-XF 

= -Lj2 dXpiX)ln[\U-iA)\] (25) 

a= —oo ^cy.q—XF 


/ OO CO 

dXp{X) [0{X — aq — Xf) 

-OO — 


-0{X-aq^XF)]\og[\UiA)\], (26) 


where 0 is the Heaviside step function and we used that 
|7n(—^)| = |/n(^)|- Furthermore, in the thermodynamic 
limit Xf = 7rn, where n is the average particle den¬ 
sity. The noncontinuous integrand will serve as the driv¬ 
ing term of the GTBA equations. Note that in the 
second line we implicitely assume that p{X) = 0 when 
A ^ [X — aq — Xf, X — aq A^?] for any a G Z. The rea¬ 
son is that for Bethe states that do not obey this condi¬ 
tion, the overlap is exactly zero (rapidities will never end 
up in those regions) and therefore 5'[p] = oo. These states 
are therefore infinitely suppressed in the Quench Action 
saddle-point equations. Another way of seeing this is 
that originally the functional integral in the quench ac¬ 
tion approach is a sum over states with non-zero overlaps 
and these states are not in that sum. 

Even when you restrict the support of the density func¬ 
tion to these intervals, in this ensemble of states there 
are still many microstates that have zero overlap with 
the Bragg-pulsed ground state. The reason is that when 
a rapidity A^^ has moved to an interval a = Pj^ it is 
not in the other intervals a ^ pj and therefore leaves a 
hole there. This alters the usual form of the Yang-Yang 
entropy significantly. Given the fillings {^a}^-cx 3 5 
finite-size entropy is 




} = 


N\ 


n“-oo(«c«!) ’ 


(27a) 


which leads to a modified Yang-Yang entropy for the 
Bragg pulse from the ground state, 

/ oo 

dXp{X)log[2Trp{X)] , (28) 

-OO 


0 = 2^ [e{X-aq-XF)-0{X-aq+XF)] log[|/„(a)|] 

a= —OO 

log(27rpgq(A)) + l-h. (29) 
This is solved by the normalized saddle-point density 


of Eq. (21), with h = 1. Eor smaller Bragg momenta 


q < 2 Xf the derivation of the saddle-point distribution 
using the QA approach remains an open problem, since 
the determinant structure of the overlaps prevents ob¬ 
taining a straightforward thermodynamic limit of the 
overlap coefficients that is expressible in terms of a root 
density p{X). 

Moreover, one could question whether the time evo¬ 
lution of simple observables obtained from the QA ap¬ 


proach using the saddle-point density of Eq. (21) is valid 


also for small Bragg momenta q < 2 Xf- However, the 
analysis of the EB mapping does not show any quali¬ 
tative differences for Bragg momenta smaller or bigger 
than 2Xf and the agreement with the time evolution of 
the QA approach is excellent for all g > 0. It therefore 
seems safe to assume that the time evolution from the 
QA approach is valid for all Bragg momenta. 


Time evolution of the momentum distribution 
function on a ring in the thermodynamic limit 


In Ref. [39] the thermodynamic limit of the matrix ele¬ 
ments of the the one-body density matrix between states 
with a countable number Uq of particle-hole excitations 
{hj Pj}]=i ^ thermodynamic state \p) was com¬ 
puted. The result is obtained by decomposing the ele¬ 
ments into a Eredholm determinant and a finite-size de¬ 
terminant accounting for the excitations. The result is 


4't(a:)4-(0) 


= + K'p) det^ [W'{hi,pj)] 

- Det(l + Rp) dfd [W{hi,pj)]\, (30) 


where p{X) is the density of rapidities of the thermody¬ 
namic state. The Eredholm determinants are denoted by 
Det, where {Kp){X^p) = K{X, p)p{p). The kernels are 
given by 


where we used the Tonks-Girardeau Bethe equation 
27r[p(A) + PhW] = 1- The variation of the quench action 
should be restricted to densities for which p{X) = 0 when 
A ^ [X — aq — Xf , X — aq Xf] for any a G Z. Also, a 


K'{X,p)=K{X,p)+ne-^^^^+>^^ , 
sin(f(A-/i)) 


K{X,p) = -An- 


X — p 


(31) 

(32) 








where n = N/L is the density. The function W is defined 
as 

W{X,fi) = {{l + Kp)-^K){X,fi), (33) 

W'{X,ii) = {{l + K'p)-^K'){X,ii), (34) 

or equivalently via the integral equations 

/ CX) 

diyK{X,v)p{iy)W{v,p) = K{X,p) , 

-OO 

(35) 

/ OO 

di/ K'{X, iy)p{iy)W'{iy, p) = K'{X, p). 

-OO 

(36) 

The QA approach yields the following expression for the 
time evolution of the one-body density matrix in the ther¬ 
modynamic limit 

V’q,A(i) 4'^(a;)4'(0) 'ipg,A{t) 

„ 1 




ne=0 


^0 K!)^ 


' '‘-e pOO 

dhjdpj 

0=1 > 


X \Pq,A 


4't(a;)4'(0) Pg,A,{hj 


(37) 


where the effective densities of holes and particles are 
given by 

, (38) 

pt\pj) = , (39) 

with the density of holes of the saddle point given by 
Pq,AiP) ~ ^ ~ Pq,A{p)- The differential overlap for 
a single particle-hole ^-^s(p)-\-6s(h) obtained by tak¬ 
ing a finite-size realization of the saddle point state 
\^q,A) -A- \pq,A) and the modified state obtained by per¬ 
forming a single particle-hole on The ratio of the 

two overlaps in the thermodynamic limit gives the differ¬ 
ential overlap 

^-Ss(p)+6s{h) ^ {i’q,A\KA,h^p) ^ 

{Wq,A\^q,A) 

The same argument gives the energy of the single 
particle-hole excitations 




(41) 


Since the overlaps only couple states such that the differ¬ 
ence of their rapidities are multiples of q, the sum over 
the excitations on the saddle point state reduces to 

pOO 

Y[L^ dhjdpj (/ii)¥’+ (Pj) 

A — -1 J—OO 


(ne!)2 

.7 = 1 




pOO 

n^/ 

i=i Pjez 


p%\hj+Pjq) 

' PgA^j + /?i«) 


(42) 


After this substitution the sum in Eq. (37) becomes the 


definition of a Fredholm determinant, and using the def¬ 
inition of the saddle-point distribution in Eq. (21), the 


expression for the time evolution of the one-body density 
matrix can be rewritten as the difference of two Fred¬ 
holm determinants of two infinite block matrices where 
each block and for any a, /3 G Z is an operator 
acting on A = [—Ai?, Ai?], 

(«'9.A(i)|«'^(a:)^'(0)|^',,^(t)) 


3? 


DetA (l<5a,/3 -|- ^a,p) 

— DetA (l<5a,/3 -|- ^oi,p)a,peZ 


(43) 


with the operators given by 


Sa,piu,v) = '^C!y\u + aq)K'{u + aq,v + {l3 + 'y)q)^^i^^^, 

7ez 


Sa,p(u,v) = '^C^\u + aq)K{u + aq,v + {j3 + 'f)q)^^p^,^. 

TGZ 

(44) 

Here u^v ^ [—^f^ and 1 is the identity operator. The 
coefficients and the function (u) are given by 


(t) _ 7 ^(iA)// 3 _[_^( —it(q^)‘^Fixqj /2 


= 


2f 




(45) 


(46) 


In order to obtain the time evolution of the momen¬ 
tum distribution h{k,t) one needs to restrict the sum 
in Eq. (37) to excitations with zero total momentum, 
namely 


i(/c, t) = 


OO ^ / Tie pOO 

n/ ‘"•^E 

ne=0^®' y=1^-00 

X (p,,a| ^'^(a:)4/(0) \pg,A, {hj -a hj + qpj]%^) 


.M,, A+\hj+pjq) 
^ PqA^j + l^jd) 


X (5 


E”ii 0: 


,o} 


(47) 


where we denoted the Fourier transform as FT{f{x)} = 
dx f{x)e~^^^. Using the identity 


f 


dv 


-iPv 


we obtain 


= (^/3,0 , 


(48) 
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Mk.t) = ft{» £ g kt. ( 14 ^,+} ■ 


«,/?(■“> '^) = X! (■“ + + a?, I' + (/3 + , 

7GZ 


7GZ 


(49) 


Time-evolved single-particle states in the trap 


Exact momentum distribution at t = 0 for arbitrary 
interactions 


The propagator for the quantum harmonic oscillator 
(Mehler kernel) is given by 


The one-body density matrix at t = 0 (after the Bragg 
pulse) is given by 


. I muj 

K{x,y;t) = X 

y 27nsm{ujt) 

—muj{x^ + y^) cos{ujt) + 2mujxy\ 

2i sin(u;t) J 

(50) 

The single particle (SP) wavefunctions after the Bragg 
pulse can then be time evolved by integrating the initial 
wavefunctions including the cosine phase with the prop¬ 
agator: 


{Ui{q,A)¥{x)i>{y)tjB{q,A)) 

= (55) 

The latter equality follows strictly from the commuta¬ 
tion relations of the Bose fields with the density and 
thus holds irrespective of interaction or geometry. For 
the case of the ring geometry, the associated momentum 
distribution function (MDF) is 

{h{k,t = 0)) = ^ d^e^’^^Io {i2Asm{q^/2)) (4'^(C)4'(0)) 

(56) 

where we defined ^ = x — y and used the integral 


/ oo 

-OO 

OO 

/3= —OO 

£(a:+|^sin(a;t))e-“«+i)7 (51) 


where i^j{x) are the groundstate harmonic eigenfunctions 

, / N 1 /mCJ\l/4 _mAU^-r-r / / - \ 

(—j e 2 Hj (£mwa;). ( 52 ) 


The result in Eq. (51) has been obtained by using the 
following two identities 


^-izcosW ^ J„(-z7e-*"'^ , (53) 

n=—OO 

j dxe-^^-y^^Hj{ax) = 

“ 54 ) 


L 

dye~^‘^^ sin(g^/2) sin{qy+q^/2) 

=/o (i2Asin(g^/2)) , (57) 

under the assumption that qLl2'K is integer. Using the 




FIG. 7. The initial MDF (t = 0) for A = 1.5, g = Stt and 
different values of the interaction strength c. The finite-c 
interactions cause a decrease of the width of the satellites but 
do not influence their relative heights. 

convolution theorem we obtain 

(n(fc, i = 0)) = ^ f{k'){h{k - 7))gs 

k' 


(58) 
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where 

f{k) = Y f dxe'^^^Io{i2Asm{qx/2)) (59) 

^ Jo 

and {h{k))GS is the MDF of the ground state. 

Using the expansion Io{z) = 
finds 


solely in the groundstate MDF {n{k ^ ql)) qs^ leading to 
a decreasing width of the peaks as one goes from the 
hard-core limit (c ^ oo) to the BEG limit (c ^ 0). In 
contrast, Eq. ( [6^ shows that their relative heights are 
completely determined by the value of A. 

Local density approximation 


(-ir 


Ioii2Asm{qx/2)) = . n 2 

n=0 ' '' 


A""sin"”(5x/2) 


= EE 

n=0 l=—n 


(-l)»+^(2n)! 
(n!)2(n — l)\{n + l)\ 


2n 

— I e 


ilqx 


(60) 


where we used the binomium to expand in plane waves. 

The order of the sums can now be interchanged. Defin¬ 
ing the coefficients 


oo 


ci{A) = Y 

n=\l\ 


(_l)n+Z(2^)! 
(n!)^(n — l)\{n + /)! 



(61) 


The local density approximation (LDA) for the gas in 
a parabolic trap amounts to replacing the value for the 
mean density in the Quench Action result for the short 
distance fluctuations with a space-dependent density pro¬ 
file corresponding to the ground state in the trap. This 
result is considerably improved when one introduces the 
classical harmonic motion of the density profile in accor¬ 


dance with the exact t = 0 MDE Eq. (63). 


In the thermodynamic limit the ground state density 
profile in a harmonic trap is given by 


Pgs(^) = (p(^))gs = — \/ mNuj — m?uj‘^x ‘^. 


(64) 


we obtain f{k) = '^iCiSk^iq. The coefficients q(A) can 
in fact be resummed and expressed in terms of a hyper¬ 
geometric function 

«(■!) = (®5^. I'l +2I‘I + ■ 

(62) 

The t = 0 post-pulse MDE can therefore be exactly ex¬ 
pressed in terms of the MDE before the pulse {n{k))GS 
as 


oo 

{n{k,t = 0))= E c;(^)(n(fc + ;g))GS • (63) 

l= — oo 

Note that this result holds for arbitrary interaction 
strength c with {h{k))GS the appropriate ground state 
MDE. The result is plotted in Eig. for different val¬ 
ues of c. The influence of the finite interactions resides 


The QA result for the time-evolved density profile yields 
PQA{x,f,n) =limth p{x) \'lpq,A{t)) 


nm 


qXpt 
cos{xqP) 


Jp{—2Asm{q‘^/3t/2m))x 

P=—oo 

sin(gAi?/3t/m) 




(65) 


with n the mean density on the ring. The result for the 
LDA in the trap then reads 

Plda{xA) = 

Ec;(^) PQA (x- — PGsix - —sin(wt)) ] 

\ ujm uom J 


( 66 ) 


where the coefficients q {A) are given in Eq. (62). 












